home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libcommon / transform.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  14.8 KB  |  559 lines
  1. /*
  2.  * transform.c
  3.  *
  4.  * Copyright (C) 1989, 1991, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * transform.c,v 4.1 1994/08/09 07:55:18 explorer Exp
  17.  *
  18.  * transform.c,v
  19.  * Revision 4.1  1994/08/09  07:55:18  explorer
  20.  * Bump version to 4.1
  21.  *
  22.  * Revision 1.1.1.1  1994/08/08  04:52:02  explorer
  23.  * Initial import.  This is a prerelease of 4.0.6enh3, or 4.1 possibly.
  24.  *
  25.  * Revision 4.0.1.1  91/09/29  15:37:06  cek
  26.  * patch1: CoordSysTransform did not detect up-Z == -1.
  27.  * 
  28.  * Revision 4.0  91/07/17  14:32:25  kolb
  29.  * Initial version.
  30.  * 
  31.  */
  32. #include "common.h"
  33.  
  34. /*
  35.  * Matrices are indexed row-first; that is:
  36.  * matrix[ROW][COLUMN]
  37.  */
  38. /*
  39.  * Allocate new structure that holds both object-to-world and
  40.  * world-to-object space transformation structures.  It probably
  41.  * should hold pointers to these structures.
  42.  */
  43. Trans *
  44. TransCreate(tr, meth)
  45. TransRef tr;
  46. TransMethods *meth;
  47. {
  48.     Trans *res;
  49.  
  50.     res = (Trans *)share_malloc(sizeof(Trans));
  51.     res->tr = tr;
  52.     res->methods = meth;
  53.     res->animated = FALSE;
  54.     res->assoc = (ExprAssoc *)NULL;
  55.     res->prev = res->next = (Trans *)NULL;
  56.     MatrixInit(&res->trans);
  57.     MatrixInit(&res->itrans);
  58.     return res;
  59. }
  60.  
  61. void
  62. TransFree(trans)
  63. Trans *trans;
  64. {
  65.     if (trans->tr)
  66.         free((voidstar)trans->tr);
  67.     free((voidstar)trans);
  68. }
  69.  
  70. void
  71. TransAssoc(trans, ptr, expr)
  72. Trans *trans;
  73. Float *ptr;
  74. Expr *expr;
  75. {
  76.     ExprAssoc *assoc;
  77.  
  78.     if (expr->timevary) {
  79.         /*
  80.          * Gotta store the sucker.
  81.          */
  82.         trans->assoc = AssocCreate(ptr, expr, trans->assoc);
  83.         trans->animated = TRUE;
  84.     } else {
  85.         *ptr = expr->value;
  86.     }
  87.     fflush(stderr);
  88. }
  89.  
  90. /*
  91.  * Allocate new transformation 'matrix'.
  92.  */
  93. RSMatrix *
  94. MatrixCreate()
  95. {
  96.     RSMatrix *res;
  97.  
  98.     res = (RSMatrix *)share_malloc(sizeof(RSMatrix));
  99.     MatrixInit(res);
  100.     return res;
  101. }
  102.  
  103. /*
  104.  * Multiply m1 and m2, copy result into "res".
  105.  */
  106. void
  107. MatrixMult(t1, t2, res)
  108. RSMatrix *t1, *t2, *res;
  109. {
  110.     register int i;
  111.     RSMatrix tmp;
  112.  
  113.     for (i = 0; i < 3; i++) {
  114.         tmp.matrix[i][0] = t1->matrix[i][0] * t2->matrix[0][0] +
  115.                      t1->matrix[i][1] * t2->matrix[1][0] +
  116.                      t1->matrix[i][2] * t2->matrix[2][0];
  117.         tmp.matrix[i][1] = t1->matrix[i][0] * t2->matrix[0][1] +
  118.                      t1->matrix[i][1] * t2->matrix[1][1] +
  119.                      t1->matrix[i][2] * t2->matrix[2][1];
  120.         tmp.matrix[i][2] = t1->matrix[i][0] * t2->matrix[0][2] +
  121.                      t1->matrix[i][1] * t2->matrix[1][2] +
  122.                      t1->matrix[i][2] * t2->matrix[2][2];
  123.     }
  124.  
  125.     tmp.translate.x = t1->translate.x * t2->matrix[0][0] +
  126.               t1->translate.y * t2->matrix[1][0] +
  127.               t1->translate.z * t2->matrix[2][0] + t2->translate.x;
  128.     tmp.translate.y = t1->translate.x * t2->matrix[0][1] +
  129.               t1->translate.y * t2->matrix[1][1] +
  130.               t1->translate.z * t2->matrix[2][1] + t2->translate.y;
  131.     tmp.translate.z = t1->translate.x * t2->matrix[0][2] +
  132.               t1->translate.y * t2->matrix[1][2] +
  133.               t1->translate.z * t2->matrix[2][2] + t2->translate.z;
  134.     MatrixCopy(&tmp, res);
  135. }
  136.  
  137. /*
  138.  * Return transformation information to map the "coordinate system"
  139.  * with the given origin, "up" vector, radius, and up axis lengths to
  140.  * one in which the "up" vector is the Z axis and the x/y/up axes
  141.  * have unit length.  This is useful for transforming a general
  142.  * form of a primitive into a canonical, Z-axis aligned, unit size
  143.  * primitive, facilitating intersection testing.
  144.  * Assumes that "up" is normalized.
  145.  */
  146. void
  147. CoordSysTransform(origin, up, r, len, trans)
  148. Vector *origin, *up;
  149. Float r, len;
  150. Trans *trans;
  151. {
  152.     RSMatrix tmp;
  153.     Vector atmp;
  154.  
  155.     ScaleMatrix(r, r, len, &trans->trans);
  156.     if (1. - fabs(up->z) < EPSILON) {
  157.         atmp.x = 1.;
  158.         atmp.y = atmp.z = 0.;
  159.     } else {
  160.         atmp.x = up->y;
  161.         atmp.y = -up->x;
  162.         atmp.z= 0.;
  163.     }
  164.     /*
  165.      * Might want to make sure that |up->z| is < 1.
  166.      */
  167.     RotationMatrix(atmp.x, atmp.y, atmp.z, -acos(up->z), &tmp);
  168.     MatrixMult(&trans->trans, &tmp, &trans->trans);
  169.     TranslationMatrix(origin->x, origin->y, origin->z, &tmp);
  170.     MatrixMult(&trans->trans, &tmp, &trans->trans);
  171.     MatrixInvert(&trans->trans, &trans->itrans);
  172. }
  173.  
  174. void
  175. TransCopy(from, into)
  176. Trans *into, *from;
  177. {
  178.     MatrixCopy(&from->trans, &into->trans);
  179.     MatrixCopy(&from->itrans, &into->itrans);
  180. }
  181.  
  182. void
  183. TransInvert(from, into)
  184. Trans *into, *from;
  185. {
  186.     RSMatrix ttmp;
  187.     /*
  188.      * In case into == from...
  189.      */
  190.     if (from == into) {
  191.         ttmp = from->trans;
  192.         into->trans = from->itrans;
  193.         into->itrans = ttmp;
  194.     } else {
  195.         into->trans = from->itrans;
  196.         into->itrans = from->trans;
  197.     }
  198. }
  199.  
  200. /*
  201.  * Copy a given transformation structure.
  202.  */
  203. void
  204. MatrixCopy(from, into)
  205. RSMatrix *into, *from;
  206. {
  207.     into->matrix[0][0] = from->matrix[0][0];
  208.     into->matrix[0][1] = from->matrix[0][1];
  209.     into->matrix[0][2] = from->matrix[0][2];
  210.     into->matrix[1][0] = from->matrix[1][0];
  211.     into->matrix[1][1] = from->matrix[1][1];
  212.     into->matrix[1][2] = from->matrix[1][2];
  213.     into->matrix[2][0] = from->matrix[2][0];
  214.     into->matrix[2][1] = from->matrix[2][1];
  215.     into->matrix[2][2] = from->matrix[2][2];
  216.     into->translate = from->translate;
  217. }
  218.  
  219. void
  220. TransInit(trans)
  221. Trans *trans;
  222. {
  223.     MatrixInit(&trans->trans);
  224.     MatrixInit(&trans->itrans);
  225. }
  226.  
  227. void
  228. TransCompose(t1, t2, res)
  229. Trans *t1, *t2, *res;
  230. {
  231.     MatrixMult(&t1->trans, &t2->trans, &res->trans);
  232.     MatrixMult(&t2->itrans, &t1->itrans, &res->itrans);
  233. }
  234.  
  235. /*
  236.  * Initialize transformation structure.
  237.  */
  238. void
  239. MatrixInit(trans)
  240. RSMatrix *trans;
  241. {
  242.     trans->matrix[0][0] = trans->matrix[1][1] = trans->matrix[2][2] = 1.;
  243.     trans->matrix[0][1] = trans->matrix[0][2] = trans->matrix[1][0] =
  244.     trans->matrix[1][2] = trans->matrix[2][0] = trans->matrix[2][1] = 0.;
  245.     trans->translate.x = trans->translate.y = trans->translate.z = 0.;
  246. }
  247.  
  248. /*
  249.  * Calculate inverse of the given transformation structure.
  250.  */
  251. void
  252. MatrixInvert(trans, inverse)
  253. RSMatrix *inverse, *trans;
  254. {
  255.     RSMatrix ttmp;
  256.     int i;
  257.     Float d;
  258.     extern int yylineno;
  259.  
  260.     ttmp.matrix[0][0] = trans->matrix[1][1]*trans->matrix[2][2] -
  261.                 trans->matrix[1][2]*trans->matrix[2][1];
  262.     ttmp.matrix[1][0] = trans->matrix[1][0]*trans->matrix[2][2] -
  263.                 trans->matrix[1][2]*trans->matrix[2][0];
  264.     ttmp.matrix[2][0] = trans->matrix[1][0]*trans->matrix[2][1] -
  265.                 trans->matrix[1][1]*trans->matrix[2][0];
  266.  
  267.     ttmp.matrix[0][1] = trans->matrix[0][1]*trans->matrix[2][2] -
  268.                 trans->matrix[0][2]*trans->matrix[2][1];
  269.     ttmp.matrix[1][1] = trans->matrix[0][0]*trans->matrix[2][2] -
  270.                 trans->matrix[0][2]*trans->matrix[2][0];
  271.     ttmp.matrix[2][1] = trans->matrix[0][0]*trans->matrix[2][1] -
  272.                 trans->matrix[0][1]*trans->matrix[2][0];
  273.  
  274.     ttmp.matrix[0][2] = trans->matrix[0][1]*trans->matrix[1][2] -
  275.                 trans->matrix[0][2]*trans->matrix[1][1];
  276.     ttmp.matrix[1][2] = trans->matrix[0][0]*trans->matrix[1][2] -
  277.                 trans->matrix[0][2]*trans->matrix[1][0];
  278.     ttmp.matrix[2][2] = trans->matrix[0][0]*trans->matrix[1][1] -
  279.                 trans->matrix[0][1]*trans->matrix[1][0];
  280.  
  281.     d = trans->matrix[0][0]*ttmp.matrix[0][0] -
  282.         trans->matrix[0][1]*ttmp.matrix[1][0] +
  283.         trans->matrix[0][2]*ttmp.matrix[2][0];
  284.  
  285.     if (fabs(d) < EPSILON*EPSILON)
  286.         RLerror(RL_PANIC, "Singular matrix.\n",yylineno);
  287.  
  288.     ttmp.matrix[0][0] /= d;
  289.     ttmp.matrix[0][2] /= d;
  290.     ttmp.matrix[1][1] /= d;
  291.     ttmp.matrix[2][0] /= d;
  292.     ttmp.matrix[2][2] /= d;
  293.  
  294.     d = -d;
  295.  
  296.     ttmp.matrix[0][1] /= d;
  297.     ttmp.matrix[1][0] /= d;
  298.     ttmp.matrix[1][2] /= d;
  299.     ttmp.matrix[2][1] /= d;
  300.  
  301.     ttmp.translate.x = -(ttmp.matrix[0][0]*trans->translate.x +
  302.                  ttmp.matrix[1][0]*trans->translate.y +
  303.                  ttmp.matrix[2][0]*trans->translate.z);
  304.     ttmp.translate.y = -(ttmp.matrix[0][1]*trans->translate.x +
  305.                  ttmp.matrix[1][1]*trans->translate.y +
  306.                  ttmp.matrix[2][1]*trans->translate.z);
  307.     ttmp.translate.z = -(ttmp.matrix[0][2]*trans->translate.x +
  308.                  ttmp.matrix[1][2]*trans->translate.y +
  309.                  ttmp.matrix[2][2]*trans->translate.z);
  310.  
  311.     MatrixCopy(&ttmp, inverse);
  312. }
  313.  
  314. /* find eigenvalues of positive definite matrices, return FALSE is singular. */
  315. int MatrixEigenvalues(trans, lambda)
  316. RSMatrix *trans;
  317. Float lambda[3];
  318. {
  319.         Float a[3], maxval, a11, a12, a21, a22;
  320.         Float *A[3];
  321.  
  322.         a[0] = fabs(trans->matrix[0][0]); maxval = a[0];
  323.         a[1] = fabs(trans->matrix[1][0]); if (a[1] > maxval) maxval = a[1];
  324.         a[2] = fabs(trans->matrix[2][0]); if (a[2] > maxval) maxval = a[2];
  325.   
  326.         if (maxval == a[0]) {
  327.                 A[0] = trans->matrix[0];
  328.                 A[1] = trans->matrix[1];
  329.                 A[2] = trans->matrix[2];
  330.         } else if (maxval == a[1]) {
  331.                 A[0] = trans->matrix[1];
  332.                 A[1] = trans->matrix[0];
  333.                 A[2] = trans->matrix[2];
  334.         } else {
  335.                 A[0] = trans->matrix[2];
  336.                 A[1] = trans->matrix[1];
  337.                 A[2] = trans->matrix[0];
  338.         }
  339.  
  340.         lambda[0] = A[0][0];
  341.         if (lambda[0] == 0.) return FALSE;
  342.  
  343.         a[1] = A[1][0] / lambda[0];
  344.         a[2] = A[2][0] / lambda[0];
  345.         a11 = A[1][1] - a[1] * A[0][1];        
  346.         a21 = A[2][1] - a[2] * A[0][1];        
  347.         a12 = A[1][2] - a[1] * A[0][2];        
  348.         a22 = A[2][2] - a[2] * A[0][2];        
  349.  
  350.         if (fabs(a11) > fabs(a21)) {
  351.                 lambda[1] = a11;
  352.                 if (lambda[1] == 0.) return FALSE;
  353.                 lambda[2] = a22 - a21/a11 * a12;
  354.         } else {
  355.                 lambda[1] = a21;
  356.                 if (lambda[1] == 0.) return FALSE;
  357.                 lambda[2] = a12 - a11/a21 * a22;
  358.         }
  359.  
  360.         if (lambda[2] == 0.) return FALSE;
  361.         return TRUE;
  362. }
  363.  
  364. /* Lipschitz number for euclidian norm: largest scaling factor */
  365. Float Lipschitz(lambda)
  366. Float lambda[3];
  367. {
  368.         Float maxlambda;
  369.  
  370.         maxlambda = lambda[0];
  371.         if (lambda[1] > maxlambda) maxlambda = lambda[1];
  372.         if (lambda[2] > maxlambda) maxlambda = lambda[2];
  373.         return maxlambda;
  374. }
  375.  
  376. /* returns TRUE if the three scaling factors are equal */
  377. int Uniform(lambda)
  378. Float lambda[3];
  379. {
  380.         return (equal(lambda[0], lambda[1]) && equal(lambda[1], lambda[2]));
  381. }
  382.  
  383. /* scaling factors */
  384. int MatrixScaling(trans, lambda)
  385. RSMatrix *trans;
  386. Float lambda[3];
  387. {
  388.         RSMatrix Q2;
  389.         int i;
  390.  
  391. /* Q-R decomposition: we are only interested in Q**2 */
  392.         Q2.matrix[0][0] = trans->matrix[0][0] * trans->matrix[0][0] +
  393.                           trans->matrix[0][1] * trans->matrix[0][1] +
  394.                           trans->matrix[0][2] * trans->matrix[0][2] ;
  395.         Q2.matrix[0][1] = trans->matrix[0][0] * trans->matrix[1][0] +
  396.                           trans->matrix[0][1] * trans->matrix[1][1] +
  397.                           trans->matrix[0][2] * trans->matrix[1][2] ;
  398.         Q2.matrix[0][2] = trans->matrix[0][0] * trans->matrix[2][0] +
  399.                           trans->matrix[0][1] * trans->matrix[2][1] +
  400.                           trans->matrix[0][2] * trans->matrix[2][2] ;
  401.         Q2.matrix[1][1] = trans->matrix[1][0] * trans->matrix[1][0] +
  402.                           trans->matrix[1][1] * trans->matrix[1][1] +
  403.                           trans->matrix[1][2] * trans->matrix[1][2] ;
  404.         Q2.matrix[1][2] = trans->matrix[1][0] * trans->matrix[2][0] +
  405.                           trans->matrix[1][1] * trans->matrix[2][1] +
  406.                           trans->matrix[1][2] * trans->matrix[2][2] ;
  407.         Q2.matrix[2][2] = trans->matrix[2][0] * trans->matrix[2][0] +
  408.                           trans->matrix[2][1] * trans->matrix[2][1] +
  409.                           trans->matrix[2][2] * trans->matrix[2][2] ;
  410.         Q2.matrix[1][0] = Q2.matrix[0][1];
  411.         Q2.matrix[2][0] = Q2.matrix[0][2];
  412.         Q2.matrix[2][1] = Q2.matrix[1][2];
  413.  
  414. /* scaling factors are square root of eigenvalues of Q**2 */
  415.         if (!MatrixEigenvalues(&Q2, lambda)) return FALSE;
  416.         lambda[0] = sqrt(lambda[0]);
  417.         lambda[1] = sqrt(lambda[1]);
  418.         lambda[2] = sqrt(lambda[2]);
  419.         return TRUE;
  420. }
  421.  
  422. /*
  423.  * Apply a transformation to a point (translation affects the point).
  424.  */
  425. void
  426. PointTransform(vec, trans)
  427. Vector *vec;
  428. RSMatrix *trans;
  429. {
  430.     Vector tmp;
  431.  
  432.     tmp.x = vec->x * trans->matrix[0][0] + vec->y * trans->matrix[1][0] +
  433.             vec->z * trans->matrix[2][0] + trans->translate.x;
  434.     tmp.y = vec->x * trans->matrix[0][1] + vec->y * trans->matrix[1][1] +
  435.             vec->z * trans->matrix[2][1] + trans->translate.y;
  436.     tmp.z = vec->x * trans->matrix[0][2] + vec->y * trans->matrix[1][2] +
  437.             vec->z * trans->matrix[2][2] + trans->translate.z;
  438.     *vec = tmp;
  439. }
  440.  
  441. /*
  442.  * 'c1x' is the X (0th) component of the first column, and so on.
  443.  */
  444. void
  445. ArbitraryMatrix(c1x, c2x, c3x, c1y, c2y, c3y, c1z, c2z, c3z, tx, ty, tz, trans)
  446. Float c1x, c1y, c1z, c2x, c2y, c2z, c3x, c3y, c3z, tx, ty, tz;
  447. RSMatrix *trans;
  448. {
  449.     trans->matrix[0][0] = c1x;
  450.     trans->matrix[1][0] = c1y;
  451.     trans->matrix[2][0] = c1z;
  452.  
  453.     trans->matrix[0][1] = c2x;
  454.     trans->matrix[1][1] = c2y;
  455.     trans->matrix[2][1] = c2z;
  456.  
  457.     trans->matrix[0][2] = c3x;
  458.     trans->matrix[1][2] = c3y;
  459.     trans->matrix[2][2] = c3z;
  460.  
  461.     trans->translate.x = tx;
  462.     trans->translate.y = ty;
  463.     trans->translate.z = tz;
  464. }
  465.  
  466. /*
  467.  * Apply transformation to a vector (translations have no effect).
  468.  */
  469. void
  470. VecTransform(vec, trans)
  471. Vector *vec;
  472. RSMatrix *trans;
  473. {
  474.     Vector tmp;
  475.  
  476.     tmp.x = vec->x*trans->matrix[0][0] +
  477.         vec->y*trans->matrix[1][0] + vec->z*trans->matrix[2][0];
  478.     tmp.y = vec->x*trans->matrix[0][1] +
  479.         vec->y*trans->matrix[1][1] + vec->z*trans->matrix[2][1];
  480.     tmp.z = vec->x*trans->matrix[0][2] +
  481.         vec->y*trans->matrix[1][2] + vec->z*trans->matrix[2][2];
  482.  
  483.     *vec = tmp;
  484. }
  485.  
  486. /*
  487.  * Transform normal -- multiply by the transpose of the given
  488.  * matrix (which is the inverse of the 'desired' transformation).
  489.  */
  490. void
  491. NormalTransform(norm, it)
  492. Vector *norm;
  493. RSMatrix *it;
  494. {
  495.     Vector onorm;
  496.  
  497.     onorm = *norm;
  498.  
  499.     norm->x = onorm.x*it->matrix[0][0] + onorm.y*it->matrix[0][1] +
  500.                 onorm.z*it->matrix[0][2];
  501.     norm->y = onorm.x*it->matrix[1][0] + onorm.y*it->matrix[1][1] +
  502.                 onorm.z*it->matrix[1][2];
  503.     norm->z = onorm.x*it->matrix[2][0] + onorm.y*it->matrix[2][1] +
  504.                 onorm.z*it->matrix[2][2];
  505.     (void)VecNormalize(norm);
  506. }
  507.  
  508. /*
  509.  * Transform "ray" by transforming the origin point and direction vector.
  510.  */
  511. Float
  512. RayTransform(ray, trans)
  513. Ray *ray;
  514. RSMatrix *trans;
  515. {
  516.         Float stretch;
  517.     PointTransform(&ray->pos, trans);
  518.     VecTransform(&ray->dir, trans);
  519.  
  520.     stretch = VecNormalize(&ray->dir);
  521.         ray->stretch *= stretch;
  522.         return stretch;
  523. }
  524.  
  525. void
  526. TransPropagate(trans)
  527. Trans *trans;
  528. {
  529.     (*trans->methods->propagate)(trans->tr, &trans->trans, &trans->itrans);
  530. }
  531.  
  532. void
  533. TransResolveAssoc(trans)
  534. Trans *trans;
  535. {
  536.     Trans *curtrans;
  537.     ExprAssoc *curassoc;
  538.  
  539.     for (curtrans = trans; curtrans; curtrans = curtrans->next) {
  540.         for (curassoc = curtrans->assoc; curassoc; curassoc = curassoc->next) {
  541.             *curassoc->lhs = ExprEval(curassoc->expr);
  542.         }
  543.         if (curtrans->assoc)
  544.             TransPropagate(curtrans);
  545.     }
  546. }
  547.  
  548. void
  549. TransComposeList(list, result)
  550. Trans *list, *result;
  551. {
  552.     TransCopy(list, result);
  553.     for (list = list->next; list; list = list->next)
  554. /* it was: 
  555.         TransCompose(list, result, result);
  556.  * PhB made it: */
  557.         TransCompose(result, list, result);
  558. }
  559.